#ifndef HITSZ_DIP_EX3_CANNY_H
#define HITSZ_DIP_EX3_CANNY_H
#include <iostream>
#include <vector>
#include <algorithm>
// #include <chrono>

#include <opencv2/opencv.hpp>
#include "MyMaxPooling2d.h"
using namespace std ;
using namespace cv;
namespace dip{
 
//******************Sobel算子计算梯度和方向********************
//第一个参数imageSourc原始灰度图像；
//第二个参数imageSobelX是X方向梯度图像；
//第三个参数imageSobelY是Y方向梯度图像；
//第四个参数pointDrection是梯度方向数组指针
//*************************************************************
void SobelGradDirction(const Mat imageSource,Mat &imageSobelX,Mat &imageSobelY,double *&pointDrection);
 
//******************计算Sobel的X和Y方向梯度幅值*************************
//第一个参数imageGradX是X方向梯度图像；
//第二个参数imageGradY是Y方向梯度图像；
//第三个参数SobelAmpXY是输出的X、Y方向梯度图像幅值
//*************************************************************
void SobelAmplitude(const Mat imageGradX,const Mat imageGradY,Mat &SobelAmpXY);
 
//******************局部极大值抑制*************************
//第一个参数imageInput输入的Sobel梯度图像；
//第二个参数imageOutPut是输出的局部极大值抑制图像；
//第三个参数pointDrection是图像上每个点的梯度方向数组指针
//*************************************************************
void LocalMaxValue(const Mat imageInput,Mat &imageOutput,double *pointDrection);
 
//******************双阈值处理*************************
//第一个参数imageInput输入和输出的的Sobel梯度幅值图像；
//第二个参数lowThreshold是低阈值
//第三个参数highThreshold是高阈值
//******************************************************
void DoubleThreshold(Mat &imageIput,double lowThreshold,double highThreshold);
 
//******************双阈值中间像素连接处理*********************
//第一个参数imageInput输入和输出的的Sobel梯度幅值图像；
//第二个参数lowThreshold是低阈值
//第三个参数highThreshold是高阈值
//*************************************************************
void DoubleThresholdLink(Mat &imageInput,double lowThreshold,double highThreshold);
 
void Canny(cv::Mat &input,cv::Mat &output,double low_thre,double high_thre)
{
    Mat imageSobelY;
	Mat imageSobelX;
    double *pointDirection=new double[(imageSobelX.cols-1)*(imageSobelX.rows-1)];  //定义梯度方向角数组
	SobelGradDirction(input,imageSobelX,imageSobelY,pointDirection);  //计算X、Y方向梯度和方向角   
    cv::Mat SobelGradAmpl;
	SobelAmplitude(imageSobelX,imageSobelY,SobelGradAmpl);   //计算X、Y方向梯度融合幅值
    Mat imageLocalMax;
	// LocalMaxValue(SobelGradAmpl,imageLocalMax,pointDirection);  //局部非极大值抑制
    maxpooling2d(SobelGradAmpl,imageLocalMax,3);
	DoubleThreshold(imageLocalMax, low_thre,high_thre);        //双阈值处理
	DoubleThresholdLink(imageLocalMax, low_thre,high_thre);   //双阈值中间阈值滤除及连接
    imageLocalMax.copyTo(output);
}
 

//******************Sobel算子计算X、Y方向梯度和梯度方向角********************
//第一个参数imageSourc原始灰度图像；
//第二个参数imageSobelX是X方向梯度图像；
//第三个参数imageSobelY是Y方向梯度图像；
//第四个参数pointDrection是梯度方向角数组指针
//*************************************************************
void SobelGradDirction(const Mat imageSource,Mat &imageSobelX,Mat &imageSobelY,double *&pointDrection)
{
	pointDrection=new double[(imageSource.rows-1)*(imageSource.cols-1)];
	for(int i=0;i<(imageSource.rows-1)*(imageSource.cols-1);i++)
	{
		pointDrection[i]=0;
	}
	imageSobelX=Mat::zeros(imageSource.size(),CV_32SC1);
	imageSobelY=Mat::zeros(imageSource.size(),CV_32SC1);
	uchar *P=imageSource.data;  
	uchar *PX=imageSobelX.data;  
	uchar *PY=imageSobelY.data;  
 
	int step=imageSource.step;  
	int stepXY=imageSobelX.step; 
	int k=0;
	int m=0;
	int n=0;
	for(int i=1;i<(imageSource.rows-1);i++)  
	{  
		for(int j=1;j<(imageSource.cols-1);j++)  
		{  
			//通过指针遍历图像上每一个像素 
			double gradY=P[(i-1)*step+j+1]+P[i*step+j+1]*2+P[(i+1)*step+j+1]-P[(i-1)*step+j-1]-P[i*step+j-1]*2-P[(i+1)*step+j-1];
			PY[i*stepXY+j*(stepXY/step)]=abs(gradY);
			double gradX=P[(i+1)*step+j-1]+P[(i+1)*step+j]*2+P[(i+1)*step+j+1]-P[(i-1)*step+j-1]-P[(i-1)*step+j]*2-P[(i-1)*step+j+1];
            PX[i*stepXY+j*(stepXY/step)]=abs(gradX);
			if(gradX==0)
			{
				gradX=0.00000000000000001;  //防止除数为0异常
			}
			pointDrection[k]=atan(gradY/gradX)*57.3;//弧度转换为度
			pointDrection[k]+=90;
			k++;
		}  
	} 
	convertScaleAbs(imageSobelX,imageSobelX);
	convertScaleAbs(imageSobelY,imageSobelY);
}
//******************计算Sobel的X和Y方向梯度幅值*************************
//第一个参数imageGradX是X方向梯度图像；
//第二个参数imageGradY是Y方向梯度图像；
//第三个参数SobelAmpXY是输出的X、Y方向梯度图像幅值
//*************************************************************
void SobelAmplitude(const Mat imageGradX,const Mat imageGradY,Mat &SobelAmpXY)
{
	SobelAmpXY=Mat::zeros(imageGradX.size(),CV_32FC1);
	for(int i=0;i<SobelAmpXY.rows;i++)
	{
		for(int j=0;j<SobelAmpXY.cols;j++)
		{
			SobelAmpXY.at<float>(i,j)=sqrt(imageGradX.at<uchar>(i,j)*imageGradX.at<uchar>(i,j)+imageGradY.at<uchar>(i,j)*imageGradY.at<uchar>(i,j));
		}
	}
	convertScaleAbs(SobelAmpXY,SobelAmpXY);
}
//******************局部极大值抑制*************************
//第一个参数imageInput输入的Sobel梯度图像；
//第二个参数imageOutPut是输出的局部极大值抑制图像；
//第三个参数pointDrection是图像上每个点的梯度方向数组指针
//*************************************************************
void LocalMaxValue(const Mat imageInput,Mat &imageOutput,double *pointDrection)
{
	//imageInput.copyTo(imageOutput);
	imageOutput=imageInput.clone();
	int k=0;
	for(int i=1;i<imageInput.rows-1;i++)
	{
		for(int j=1;j<imageInput.cols-1;j++)
		{
			int value00=imageInput.at<uchar>((i-1),j-1);
			int value01=imageInput.at<uchar>((i-1),j);
			int value02=imageInput.at<uchar>((i-1),j+1);
			int value10=imageInput.at<uchar>((i),j-1);
			int value11=imageInput.at<uchar>((i),j);
			int value12=imageInput.at<uchar>((i),j+1);
			int value20=imageInput.at<uchar>((i+1),j-1);
			int value21=imageInput.at<uchar>((i+1),j);
			int value22=imageInput.at<uchar>((i+1),j+1);
 
			if(pointDrection[k]>0&&pointDrection[k]<=45)
			{
				if(value11<=(value12+(value02-value12)*tan(pointDrection[i*imageOutput.rows+j]))||(value11<=(value10+(value20-value10)*tan(pointDrection[i*imageOutput.rows+j]))))
				{
					imageOutput.at<uchar>(i,j)=0;
				}
			}	
			if(pointDrection[k]>45&&pointDrection[k]<=90)
 
			{
				if(value11<=(value01+(value02-value01)/tan(pointDrection[i*imageOutput.cols+j]))||value11<=(value21+(value20-value21)/tan(pointDrection[i*imageOutput.cols+j])))
				{
					imageOutput.at<uchar>(i,j)=0;
 
				}
			}
			if(pointDrection[k]>90&&pointDrection[k]<=135)
			{
				if(value11<=(value01+(value00-value01)/tan(180-pointDrection[i*imageOutput.cols+j]))||value11<=(value21+(value22-value21)/tan(180-pointDrection[i*imageOutput.cols+j])))
				{
					imageOutput.at<uchar>(i,j)=0;
				}
			}
			if(pointDrection[k]>135&&pointDrection[k]<=180)
			{
				if(value11<=(value10+(value00-value10)*tan(180-pointDrection[i*imageOutput.cols+j]))||value11<=(value12+(value22-value11)*tan(180-pointDrection[i*imageOutput.cols+j])))
				{
					imageOutput.at<uchar>(i,j)=0;
				}
			}
			k++;
		}
	}
}
 
//******************双阈值处理*************************
//第一个参数imageInput输入和输出的的Sobel梯度幅值图像；
//第二个参数lowThreshold是低阈值
//第三个参数highThreshold是高阈值
//******************************************************
void DoubleThreshold(Mat &imageIput,double lowThreshold,double highThreshold)
{
	for(int i=0;i<imageIput.rows;i++)
	{
		for(int j=0;j<imageIput.cols;j++)
		{
			if(imageIput.at<uchar>(i,j)>highThreshold)
			{
				imageIput.at<uchar>(i,j)=255;
			}	
			if(imageIput.at<uchar>(i,j)<lowThreshold)
			{
				imageIput.at<uchar>(i,j)=0;
			}	
		}
	}
}
//******************双阈值中间像素连接处理*********************
//第一个参数imageInput输入和输出的的Sobel梯度幅值图像；
//第二个参数lowThreshold是低阈值
//第三个参数highThreshold是高阈值
//*************************************************************
void DoubleThresholdLink(Mat &imageInput,double lowThreshold,double highThreshold)
{
	for(int i=1;i<imageInput.rows-1;i++)
	{
		for(int j=1;j<imageInput.cols-1;j++)
		{
			if(imageInput.at<uchar>(i,j)>lowThreshold&&imageInput.at<uchar>(i,j)<255)
			{
				if(imageInput.at<uchar>(i-1,j-1)==255||imageInput.at<uchar>(i-1,j)==255||imageInput.at<uchar>(i-1,j+1)==255||
					imageInput.at<uchar>(i,j-1)==255||imageInput.at<uchar>(i,j)==255||imageInput.at<uchar>(i,j+1)==255||
					imageInput.at<uchar>(i+1,j-1)==255||imageInput.at<uchar>(i+1,j)==255||imageInput.at<uchar>(i+1,j+1)==255)
				{
					imageInput.at<uchar>(i,j)=255;
					DoubleThresholdLink(imageInput,lowThreshold,highThreshold); //递归调用
				}	
				else
			{
					imageInput.at<uchar>(i,j)=0;
			}				
			}				
		}
	}
}
 
}

#endif 